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ABSTRACT 

The direct simulation Monte-Carlo (DSMC) method 
was applied to the analysis of low-density nitrogen plumes 
exhausting from a small converging-diverging nozzle into 
finite ambient pressures. Two cases were considered that 
simulated actual test conditions in a vacuum facility. The 
numerical simulations readily captured the complicated 
flow structure of the overexpanded plumes adjusting to the 
finite ambient pressures, including Mach disks and barrel- 
shaped shocks. The numerical simulations compared well 
to experimental data of Rothe. 

INTRODUCTION 

Low-thrust rocket engines which are used for station- 
keeping, attitude and altitude control on various satellites 
and spacecraft have a significant impact on mission 
performance such as on-orbit lifetime, payload, and trip 
time. Another important factor affecting mission 
performance is contamination of sensitive instruments and 
system components that may be exposed to the thruster 
plumes both in the forward and backward flow directions. 
Hence, understanding of the detailed flow structure in and 
around low-thrust rocket nozzles is important not only for 
the accurate prediction of thrust and mass flow but also for 
the precise prediction of the plume flow regions. 

Low-density flow through small nozzles expanding to 
vacuum has been examined previously in both 
experimental and numerical investigations, though, little 
experimental data are available and most data are for gross 
characteristics of nozzle performance such as thrust and 
discharge coefficients. Data that provides detailed 
information on internal flow structure was published by 
Rothe [1,2] in which density and rotational temperature 
distributions were measured inside a small nozzle and in the 
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nozzle plume using the electron-beam fluorescence 
technique. Rothe's experiment was numerically simulated 
by Chung et al [3] using a continuum code based on the 
Navier-Stokes equations and the direct simulation Monte- 
Carlo (DSMC) method of Bird [4], The simulation results 
were compared with Rothe's density and rotational 
temperature data at various locations inside the nozzle and 
at the nozzle exit plane. In this study, no consideration 
was given to simulation of the finite test-facility pressure. 
In somewhat related work, Pitot-pressure measurements and 
numerical simulations were made for a low-density nozzle 
and plume flow by Penko et al [5] and Boyd et al [6]. 
Comparisons were made between continuum and DSMC 
results and Pitot pressure measurements at the nozzle exit 
plane and various locations in the plume. Comparisons of 
continuum and DSMC results were also made for the flow 
inside the nozzle. In these works, test-facility pressures 
were quite low and, therefore, consideration was not given 
to numerical simulation of the actual ambient pressure. In 
other work, Campbell [7,8], Nelson and Doo [9], and 
Zelesnik et al [10] also analyzed expanding low-density 
nozzle flows using the DSMC method and compared their 
results with experimental data. 

The rocket engines under consideration typically have 
thrust levels under 100 mN, are physically small in size, 
and have relatively low operating pressures and mass flow. 
Reynolds numbers are low, on the order 10 3 , and 
rarefaction effects are significant. Under these conditions, 
the flow contains strong nonequilibrium effects, such as 
slip at a wall, from rapid expansion to near vacuum 
conditions. The flow transitions from continuum to free- 
molecular regimes. Conventional continuum gas dynamics 
may not be adequate to accurately analyze the flow and an 
approach based on kinetic theory may be required. 

Of the various methods available for analysis of low- 
density gas flows, the DSMC method is most widely used 
and readily applicable. The DSMC method is a numerical 
simulation technique for solving the Boltzmann equation 
by modeling a real gas flow using a representative set of 
molecules. Theoretically, the DSMC method can be 
applied to any flow for which the Boltzmann equation is 
valid but intensive computational requirements generally 
restrict the use to near continuum and rarefied flows. 
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Continuum methods are usually much more efficient than 
the DSMC method for higher density flows. Thus, in the 
analysis of flows which involve both continuum and 
rarefied regimes, it is reasonable to apply both methods. 
The simplest utilization of both methods is to solve the 
rarefied flow regime with the DSMC method using 
boundary conditions for the inflow surface obtained from 
the continuum method used to solve the continuum and 
near-continuum regime. 

In this paper, the DSMC method is employed for the 
analysis of low-density nitrogen flows expanding through a 
small converging-diverging nozzle and into finite back 
pressures that simulate Rothe's experimental conditions. 
Special attention is given to the effect of the non-zero 
ambient pressure on the flow structure, both in the nozzle 
but especially in the plume, to simulate conditions often 
encountered in ground-based test facilities used to test 
small thrusters. In contrast to the behavior of low-density 
plumes expanding into a vacuum, the flow expanding into a 
region of finite pressure is often overexpanded with a more 
complicated flow pattern involving Mach disks and shock 
waves. Under these conditions, the supersonic flow can be 
confined to a narrow core, depending on the applied 
background pressure, and the sonic line may not intersect 
with the nozzle lip as it does with a vacuum ambient, in 
which case the external conditions can influence the 
internal flow through the thick subsonic region near the 
nozzle wall. These conditions are illustrated in the results 
from the numerical simulation and in comparison with 
Rothe's experimental data. 

PROBLEM STATEMENT 

Rothe’s experiment [2] was chosen as a reference 
problem because of the availability of detailed 
measurements. Figure 1 illustrates the geometry of the 
nozzle used in Rothe's experiment and in the numerical 
simulation. The actual nozzle was made of graphite to 
reduce optical reflections and to minimize back-scattering 
and secondary emission of electrons. The subsonic and 
supersonic portions of the nozzle are cones having half- 
angles of 30° and 20°, respectively, with longitudinal radii 
of curvature at the throat equal to 1/2 of the throat radius. 
The area ratio at the exit based on the throat area is 66. 
The shaded region in Fig. 1 indicates the domain of the 
DSMC simulation. The length of the curved contour 
downstream of the throat (IH) is about 0.5 mm. The 
simulation domain extends from the throat to an axial 
distance 260 mm from the nozzle exit plane (AB), a radial 
distance 50 mm from the nozzle lip (GD), and an axial 
distance 20 mm from the nozzle exit plane into the 
backflow region (DE). The inflow boundary is located at 
the nozzle throat (OI). The boundary condition at the 
inflow boundary is obtained from a solution of the Navier- 
Stokes equations [3]. The radial boundary CE is located 
far enough from the axis so that extending it further 
radially does not result in any significant changes in the 
macroscopic flow variables in the plume. Along the 
boundaries CE and EF an equilibrium condition 
corresponding to the ambient gas is used as the boundary 


condition. Likewise, the downstream boundary BC is 
located far enough from the nozzle exit so that extending it 
further downstream does not result in any significant 
changes in the macroscopic flow variables in the near 
plume. Along the boundary BC, several boundary 
conditions are tested including an equilibrium condition 
corresponding to the ambient condition, a vacuum 
boundary condition, and an equilibrium condition 
corresponding to a profile extrapolated from the inside of 
the plume. The test gas is nitrogen with a stagnation 
temperature of T 0 = 300 K. The flow conditions are listed 
in Table 1. In the table, the throat Reynolds number, 
fl e ,t = 2ml%\i 0 Ru is based on the viscosity at the stagnation 
chamber condition, |i 0 . Here the quantity m is the mass 
flow rate and Rt is the throat radius. The Knudsen number 
Kn is based on the throat diameter and the stagnation 
chamber condition. 

DSMC METHOD 

The DSMC code used in this study is based on the 
principles described by Bird [4], together with the variable 
hard sphere (VHS) model [11] as a molecular model and 
the no time counter (NTC) method [12] as a collision 
sampling technique. The code was developed to 
investigate various low-density flows of gas mixtures in 
arbitrarily shaped flow domains [3,13,14]. Details of the 
code may be found in Ref. 3. The execution speed of the 
code for the flow considered in this study, measured by 
CPU time/particle/timestep, is about 1.3|is on a CRAY 
Y/MP. The flow domain consists of about 20,000 cells in 
41 subregions. At the steady phase of the simulation, the 
total number of simulated molecules in the flow domain is 
about 2 million. The flow field is sampled every 5 
timesteps during 20,000 timesteps after reaching the steady 
phase. The total CPU time required for the computation is 
about 19 hours on the CRAY Y/MP. 

The YHS exponent, to, of nitrogen is chosen to be 0.24 
with the reference molecular diameter of 4.07xl(T 10 m at 
the reference temperature 273K [11]. Chemical reactions 
and the vibrational mode are assumed to be frozen. For the 
calculation of rotational energy exchange between the 
colliding molecules, the Borgnakke-Larsen phenomeno- 
logical model [15] is employed together with the 
temperature-dependent energy exchange probability of 
Boyd [16] modified by Chung et al [3] to be consistent 
with the experimental data for the rotational relaxation of 
nitrogen obtained by various methods and compatible with 
the VHS model. A diffusely adiabatic wall with 10% 
thermal accommodation is assumed for the interaction 
between the gas molecules and the wall [3]. 

To assess the effect of backscattering of downstream 
molecules and reduce the number of simulated molecules, 
the gas is treated as a mixture consisting of two different 
nitrogen sources whose origins are the upstream reservoir 
and the downstream, respectively. A ratio of the actual to 
the simulated molecules for the nitrogen whose origin is 
the downstream is chosen to be 4 times larger than that for 
the nitrogen coming from the upstream reservoir. 
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RESULTS AND DISCUSSION 

To present the general idea of the over-expanded 
plumes and the overall structure of the flowfield, density 
and Mach number contours will be considered first. 
Figures 2 and 3 show density and Mach number contours 
for Case I in which the stagnation pressure Po is 474 Pa and 
the ratio of stagnation to ambient pressure is PJPb — 310. 
The density is normalized by the background density pb. In 
this case, the density decreases from the throat to the 
nozzle exit and approaches the background density in the 
plume. Around the nozzle lip, the density has a minimum 
value which is about 20% lower than the background 
density due to entrainment. The Mach number increases 
from the throat to the nozzle exit and then decreases as the 
flow adjusts to the ambient pressure. The maximum Mach 
number is about 3.7 and occurs just inside the nozzle near 
the exit. In contrast to the behavior observed in plumes 
expanding into a vacuum, the sonic line of the internal 
boundary layer does not intersect with the nozzle lip and is 
almost parallel to the nozzle surface. There is no sign of 
shock waves in this plume. 

Figures 4 and 5 show density and Mach number 
contours for Case II in which the stagnation pressure is 
1245 Pa and pressure ratio is PJPb= 325. Again, the 
density is normalized by the background density pb. In this 
case, the plume is overexpanded and a more complicated 
flow pattern exists with a barrel-shaped shock wave and a 
Mach disk. It can be seen that the supersonic portion of the 
flow is restricted to a narrow core and enveloped by the 
barrel-shaped shock wave. The flow in the core region is 
strongly decelerated at the Mach disk and there is a region 
of high density between the barrel-shaped shock and the 
plume boundary. Also, the plume shows a repetitive 
pattern of expansion, which is diffused far downstream. 

The pressure distribution along the axis of the nozzle is 
depicted in Fig. 6, which shows the extent of 
overexpansion in the plume. It can be seen in the figure 
that even in Case I the flow is slightly overexpanded just 
outside the nozzle exit, where the minimum pressure is 
about 73% of the ambient pressure. However, the 
overexpansion is not strong enough to cause any significant 
disturbances in the flow. In case II, the centerline pressure 
at the exit is already 60% of the ambient pressure and drops 
to 33% of the ambient pressure at 25 mm downstream 
from the exit. Hence, back-pressure effects cause the flow 
to turn into the axis and produce the characteristic barrel- 
shock configuration. Beyond the shock, the flow continues 
to expand. In this case, however, the expansion is not 
sufficiently strong to sustain the repetitive pattern of 
compression and expansion. 

Consideration is now given to the effect of boundary 
conditions. For the simulation of nozzle flows with non- 
zero ambient pressures a supersonic outflow boundary 
condition, which is usually assumed along boundaries in 
the simulation of nozzle flows expanding into a vacuum, 
can not be used since the flow is not supersonic along the 
boundaries in the plume. The best way would be to extend 
the simulation domain sufficiently far so that the flow 
along the boundaries is not disturbed by the expanding 
plume and an equilibrium condition corresponding to the 


ambient gas may be used as the boundary condition, which, 
unfortunately, is computationally prohibitive. In the 
present study, the extension of the simulation domain to a 
radial distance 70 mm from the axis, and to an axial 
distance 20 mm from the nozzle exit plane into the 
backflow region, is found to be sufficient to impose the 
equilibrium condition corresponding to the far-field 
ambient condition. It is also found that the extension of the 
simulation domain to the downstream axial distance 260 
mm from the nozzle exit plane is sufficient to resolve 
important characteristics of the plume. 

Figure 7 shows the effect of the boundary condition 
along the downstream boundary on the centerline Mach 
number. Two distinctively different boundary conditions 
are imposed along the downstream boundary at z = 260 
mm, including, an equilibrium condition corresponding to 
the ambient condition (dot-dashed line) and a vacuum 
boundary condition (dashed line). For a third condition, 
since the radial distribution of flow variables are very 
similar far downstream, profiles are extrapolated from the 
inside of the plume and an equilibrium condition 
corresponding to the profiles is employed as the boundary 
condition (solid line). For the sake of comparison, the 
result obtained by reducing the axial distance from the 
nozzle exit to the downstream boundary to about 180 mm 
and imposing an equilibrium condition corresponding to 
the ambient condition is also shown (dotted line). It can be 
seen that extremely different boundary conditions have no 
significant effect on the Mach number distribution in the 
near plume within 150 mm from the nozzle exit. Figure 8 
shows the effect of the boundary condition on the radial 
translational temperature profile at the axial location z = 
150 mm. Again it can be seen that extremely different 
boundary conditions along the downstream boundary have 
no significant effect on the radial temperature distribution 
in the near plume within 150 mm from the nozzle exit, 
where the most important and distinctive characteristics of 
the present overexpanded plume occurs. Although not 
shown, the effect of the boundary conditions on the other 
flow variables such as density, velocity, and rotational 
temperature is similar to those observed in Figs. 7 and 8. 

The comparison of the simulation results with Rothe s 
experimental data[2] is made in Fig. 9 in which density 
profiles along the nozzle axis inside the nozzle and in the 
plume for Case I are shown. The densities are normalized 
by the ambient density pt>. The solid line is from the 
DSMC simulation with the ambient pressure and the filled 
circles represent the experimental data of Rothe. For 
comparison, a calculation made for the flow expanding into 
a vacuum is also shown in the figure and is represented by 
the dotted line. The simulation results indicate that for the 
two different ambient conditions, the density profiles along 
the axis show no significant differences inside the nozzle 
for the case where the stagnation pressure P 0 = 474 Pa. 
The simulation result with the non-zero ambient pressure 
shows good agreement with the experimental data both 
inside the nozzle and in the plume. Figure 10 shows the 
comparison of the simulation results with experimental 
data along the nozzle axis for Case II where P 0 = 1245 Pa. 
The density decreases from the throat to the nozzle exit 
before the shock occurs, the local minimum density is 
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about 2.2 times the ambient density and occurs about 25 
mm downstream from the exit. The local maximum 
density is about 6.4 times the ambient density and occurs 
about 60 mm downstream from the exit. The density 
profile along the axis with a non-zero ambient pressure and 
that for a vacuum ambient show no significant differences 
before the shock for this case. The simulation result with 
the non-zero ambient pressure shows good agreement with 
the experimental data across the shock except for the slight 
shift in the axial direction. 

Comparison of calculated radial density profiles with 
Rothe's experimental data at various axial locations in the 
plume for Cases I and II are shown in Figs. 1 1 and 12, 
respectively. Again, the solid line is the simulation result 
with non-zero ambient pressure and the filled circles 
represent the experimental data. The densities are 
normalized by the ambient density pb. The radial density 
profiles are arranged in a sequence corresponding to axial 
positions in the plume. At the lower stagnation pressure of 
p 0 = 474 Pa, the density profiles have a single peak, which 
decrease both in the radial and axial directions. At the 
higher stagnation pressure of Po— 1245 Pa, the density 
profiles in the plume are double peaked and indicate the 
presence of a barrel-shaped shock. The barrel-shaped 
shock converges to a single density peak further 
downstream. In both cases, the jet is confined to a narrow 
core the diameter of which is about the same as the nozzle 
exit diameter. The simulation results show good agreement 
with the experimental data in both cases except for the 
region around the shock in the higher pressure case where a 
small difference in the axial location results in a significant 
difference in the profile as shown in Fig. 10. 

The degree of thermal nonequilibrium in the plume is 
shown in Fig. 13 where both translational and rotational 
temperature distributions along the axis for Case II are 
plotted. For comparison, temperature distributions for the 
flow expanding into a vacuum are also shown in the figure. 
From the rapid expansion, the rotational temperature is 
higher than the translational temperature along the axis 
inside the nozzle and in the near plume. In the case of the 
flow expanding into a vacuum, the two temperatures are 
frozen far downstream. For the case of the non-zero 
ambient pressure, however, the translational temperature 
becomes higher than the rotational temperature across the 
shock wave from the rapid compression. Beyond the 
shock, the translational temperature is higher than the 
rotational temperature where the flow once again expands. 
Far downstream, the two temperatures become the same as 
the flow approaches equilibrium. 

Finally, the effect of the non-zero ambient pressure on 
the exit plane density is shown in Figs. 14 and 15. In the 
figures, the quantity Re is the exit radius of the nozzle. 
Figure 14 shows the exit plane density profiles for the 
stagnation pressure P 0 = 1245 Pa. Near the axis, the two 
densities are nearly the same and are not influenced by the 
external condition, while near the nozzle lip, the effect of 
the external condition becomes significant. At the nozzle 
wall, the density of the flow expanding into vacuum is 55% 
of that with the non-zero ambient pressure. Figure 15 
shows the mole fraction of the background gas at the exit 
plane. It can be seen that near the axis the gas is composed 


of gas molecules exclusively coming from the upstream 
stagnation chamber. Near the nozzle lip, a significant 
number of gas molecules from the ambient have penetrated 
inside the nozzle. At the lip, the fraction of gas originating 
from the background gas is about 11% and 13% for Case I 
and II, respectively. 

CONCLUSIONS 

The flow structure in low-density plumes expanding 
through a small converging-diverging nozzle into a region 
of finite pressure is investigated using the direct simulation 
Monte-Carlo (DSMC) method based on molecular gas 
dynamics. This kind of flow structure has rarely been 
studied using the DSMC method. The calculated results 
show that the plume structure with a non-zero ambient 
pressure is substantially different from that of the flow 
expanding into a vacuum. The simulation results show 
good agreement with Rothe's experimental data measured 
at various axial and radial locations in the plume. It is 
shown that the DSMC method can successfully predict the 
complex flow structures such as a Mach disk and barrel - 
shaped shock which are representative of overexpanded 
nozzle plumes. 
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Fig. 1 Geometry of low-thrust nozzle. 


Table 1 Flow conditions 



Case I 

Casell 

Test gas 

n 2 

n 2 

Stagnation temperature, T 0 

300 K 

300 K 

Stagnation pressure, P 0 

474 Pa 

1245 Pa 

Ambient pressure, Pb 

1.5 Pa 

3.8 Pa 

Wall temperature, 7w 

300 K 

300 K 

Reynolds number, Re , t 

270 

709 

Knudsen number, Kn 

2.3xl0" 3 

8.8xl0' 4 
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Fig. 3 Mach number contours for P 0 = 474 Pa and PJP b - 310. 



Fig. 4 Density contours for P 0 = 1245 Pa and PJP\> - 325. 
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Fig, 5 Mach number contours for P 0 = 1245 Pa and PJP h = 325. 



Fig. 6 Pressure distribution along the axis. 



Fig. 8 Translational temperature distribution at z = 150 mm. 



Fig. 7 Mach number distribution along the axis. 



Fig. 9 Density distribution along the axis for P 0 = 474 Pa. 


Normalized density, p/p. ** era Normalized density, p/p. era* Normalized density, p/p | 



Distance from nozzle exit, mm Distance from axis, mm 

10 Density distribution along the axis for P 0 = 1245 Pa, Fig. 11 Comparison of radial density distribution for 

Pq = 474 Pa. 



Distance from axis, mm 
12 Comparison of radial density distribution for 
1245 Pa. 



Normalized distance from axis, r/R £ 


Fig. 14 Density distribution at exit plane for P Q = 1245 Pa. 



Distance from nozzle exit, mm 
Fig. 13 Temperature distribution along the axis. 



Normalized distance from axis, r/R^ 
Fig. 15 Mole fraction of ambient gas at exit plane. 
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